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1.0  SUMMARY 


In  this  report,  a  discontinuous  Galerkin  time-domain  (DGTD)  method  is  developed  to  simulate 
high-power  microwave  (HPM)  air  breakdown  phenomena,  which  are  modeled  by  a  coupled 
electromagnetic— plasma  system.  In  the  coupled  system,  the  electromagnetic  fields  are  governed 
by  Maxwell’s  equations  and  the  plasma  is  modeled  by  the  five-moment  fluid  equations  (Euler’s 
equations).  The  non-Maxwellian  electron  energy  distribution  function  (EEDF)  is  used  to 
calculate  electron  transport  coefficients  and  describe  the  non-equilibrium  collision  reactions 
between  electrons  and  neutral  air  particles.  The  coupled  Maxwell— Euler  equations  are  solved  by 
the  DGTD  method  with  high-order  spatial  and  temporal  discretizations,  which  are  able  to 
provide  a  sufficient  resolution  for  the  physical  quantities  in  both  space  and  time.  Several 
numerical  examples  are  presented  to  investigate  the  physical  process  and  demonstrate  the 
capability  of  the  numerical  method. 

2.0  INTRODUCTION 

When  HPM  pulses  travel  through  air,  the  high  energy  of  the  pulses  can  ionize  neutral  air 
particles  and  cause  air  to  break  down.  In  this  process,  a  huge  number  of  free  electrons  are  generated 
and  heated  up  by  the  power  supplied  by  the  electromagnetic  fields.  The  electron  oscillation 
generates  secondary  fields  which,  when  strong  enough,  can  cancel  a  part  of  the  HPM  pulses  being 
transmitted  and  cause  the  so-called  HPM  tail  erosion  and  pulse  shortening.  This  will  severely  limit 
the  transmission  of  the  HPM  pulses  in  air  and  has  been  investigated  by  many  researchers  [1]-[1 1], 

To  describe  the  HPM  breakdown  process,  different  models  can  be  employed.  Based  on  the 
fact  of  the  insulator-to-conductor  transition  during  breakdown,  a  nonlinear  conductivity  model 
was  developed  in  [12],  where  the  material  conductivity  was  assumed  to  be  a  nonlinear  function 
of  the  electric  field.  In  [13],  a  simplified  plasma  model  considering  particle  collisions  was 
employed.  In  this  model,  the  plasma  velocity  was  treated  as  an  unknown  quantity  governed  by 
the  momentum  transfer  equation  with  collision.  The  resulting  coupled  system  was  solved  by  a 
nonlinear  finite-element  time-domain  (FETD)  method.  To  include  more  physics,  a  plasma 
diffusion  model  was  adopted  and  solved  by  a  coupled  DGTD  method  in  [14],  where  the  air 
breakdown  and  tail  erosion  were  simulated  by  considering  both  the  electron  density  and  velocity 
as  unknown  quantities.  To  obtain  a  numerical  solution  with  a  higher  fidelity,  a  five-moment 
plasma  fluid  model  is  adopted  in  this  work,  where  the  plasma  density,  velocity,  and  energy  are 
all  considered  and  described  by  three  conservation  laws  obtained  by  taking  the  first  three 
moments  of  the  Boltzmann  equation.  To  describe  the  non-equilibrium  collision  process  between 
electrons  and  neutral  particles,  the  transport  coefficients  are  calculated  by  integrating  the 
collision  cross  sections  with  the  non-Maxwellian  EEDF  [9],  [15], 

In  the  simulation  of  the  coupled  electromagnetic— plasma  system,  the  finite-difference  time- 
domain  (FDTD)  method  [1 1],  [16]  is  most  widely  used  because  of  its  simplicity  and  high  parallel 
efficiency.  However,  the  stair-case  approximation  of  the  solution  domain,  the  finite  difference 
approximation  of  the  fields  and  their  derivatives,  and  the  explicit  leap-frog  time-marching  scheme 
used  in  the  FDTD  method  result  in  an  overall  low-order  accuracy,  which  requires  an  extremely 
dense  mesh  grid  and  an  extremely  tiny  time  step  size  in  a  simulation  to  achieve  a  desired  accuracy. 
To  overcome  these  issues,  a  coupled  high-order  DGTD  scheme  [17]-[21]  is  adopted  in  this  report 

1 


Approved  for  public  release;  distribution  is  unlimited. 


to  solve  the  coupled  system  equations.  Using  an  unstructured  mesh,  high-order  nodal  basis 
functions  and  an  explicit  high-order  time-marching  scheme,  the  employed  DGTD  method  is  able 
to  achieve  a  high-order  accuracy  in  the  spatial  and  temporal  representations  of  the  physical 
quantities.  Also,  it  is  able  to  preserve  the  continuity  of  all  components  of  the  electromagnetic 
fields  due  to  the  application  of  the  nodal  basis  functions  [22],  which  is  critical  in  maintaining  the 
stability  of  a  Boltzmann  solver.  Moreover,  the  flexibility  and  efficiency  of  the  DGTD  method 
can  be  enhanced  through  the  application  of  the  dynamic  h-  and/or  /^-adaptation  techniques  [23], 
[24]  and  massive  parallelization  on  a  multi-core  CPU  cluster  or  a  many-core  graphic  processing 
unit  (GPU)  platform  [25],  [26], 

This  report  is  organized  as  follows.  In  Section  3.0,  the  five-moment  fluid  model  is  first 
introduced,  followed  by  the  description  of  the  calculation  of  its  source  terms.  The  DGTD 
solution  of  the  coupled  Maxwell— Euler  system  is  then  presented.  The  electromagnetic— plasma 
interactions  and  the  resulting  air  breakdown,  electromagnetic  pulse  tail  erosion,  plasma 
formation  and  shielding  are  simulated  and  given  as  numerical  examples  in  Section  4.0  to 
demonstrate  the  physical  process  and  validate  the  numerical  method.  Conclusions  are  drawn  in 
Section  5.0.  Note  that  in  this  report,  only  electrons  are  considered.  The  positive  and  negative 
ions  are  considered  as  immobile  due  to  their  large  mass  and  therefore,  the  ion  currents  and  their 
radiation  are  neglected. 

3.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

3.1  Plasma  Model 

3.1.1  The  Boltzmann  &  General  Transport  Equations 

The  behaviors  of  electrons  are  governed  by  the  Boltzmann  equation  [27]-[30],  which  reads 

df  F  (df\ 

-Ui '-V/  +  — -Vv/=  -M  (1) 

Ot  TTle  Vc/t/coll 

where  /  =  f(r,  v,  t)  denotes  the  electron  distribution  function  defined  in  the  seven-dimensional 
phase  space  (r,  v,  t),  me  denotes  the  electron  mass  at  rest,  and  F  denotes  the  macroscopic  force 
such  as  the  electromagnetic  (Lorentz)  and  gravitational  forces.  In  the  rest  of  this  report,  F 
denotes  the  Lorentz  force  only,  since  the  gravitational  force  is  negligible  when  comparing  with 
the  Lorentz  force  in  problems  considered  here.  The  operators  V  and  Vv  denote  the  gradients 
taken  in  the  physical  space  r  =  (x,  y,  z)  and  the  velocity  space  v  =  (  v  x,  v  y,  vz),  respectively. 
The  right-hand  side  (RHS)  of  the  Boltzmann  equation  is  an  abstract  form  representing  the 
temporal  variation  of  the  electron  distribution  function  due  to  collisions. 

In  fluid  models,  the  behaviors  of  various  discharge  particle  species  are  described  in  terms  of 
average,  macroscopic,  hydrodynamic  quantities  such  as  the  particle  density 

n(r,  t)  =  [  f(r,  v,  t)  dv  (2) 


the  mean  velocity 


2 


Approved  for  public  release;  distribution  is  unlimited. 


U(r,t )  =  (v) 


=  ~  [  vf  (r,  v, 

nj 


t)  dv 


(3) 


and  the  mean  energy 


1  771  f 

£(r,t)  =  -me(v2)  =  ^  I  v2f(r,v,t)  dv 


(4) 


In  general,  a  macroscopically  averaged  quantity  can  be  defined  as 

/  <t>(v)  f(r,v,  t )  dv 


1  f 

(d>(v))  =  —  d>(v)  f(v,  v>  t)  dv  = 

nj 


/  f(r,v,t )  dv 


(5) 


where  <t>  is  some  function  of  velocity,  which  can  be  a  scalar,  a  vector,  or  a  tensor.  In  (3)  and  (4), 
d>  are  taken  as  the  microscopic  velocity  v  and  the  microscopic  energy  ~mev2  of  a  single  electron, 
respectively. 

Multiplying  the  Boltzmann  equation  by  <J>(v)  and  integrating  over  all  velocity  components 
yield  the  general  transport  equation 


[  d>(v)  dv  +  [  <P(v)  v  ■  V/  dv  +  [  <t>(v) - Vvf  dv  =  [  d>(v)  (r£-\  dv 

J  dt  J  J  77ie  J  \ot  /  CqH 


(6) 


By  invoking  the  definition  of  the  average  quantity  (5),  the  first  term  on  the  left-hand  side  (LHS) 
of  (6)  can  be  rewritten  as 


C  d  f  d  f  d 

J  <f>  dv  =  —J  (pfdv  =  —  [n<<!>0)>] 

Using  the  integration  by  parts,  the  second  term  of  (6)  is  expressed  as 

J  d>  v  ■  Vf  dv  =  J  [V  ■  (fPvf)  -  f  V  ■  (d>v)]  dv 

Since  d>v  is  not  a  function  of  r,  V  ■  (<E>v)  =  0.  The  second  term  becomes 

J  d>  v  ■  Vf  dv  =  J  V  ■  (d>vf)  dv  =  V  ■  J  (Pvf  dv  =  V  ■  [n(d>(r)v)] 


(7) 


(8) 


(9) 


Applying  integration  by  parts  and  the  Gauss  divergence  theorem,  the  third  term  of  (6)  can  be 
rewritten  as 
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/  kVJdv  =  i  K  ■  (*£')  (*5] 

=  -{<I>-/d5p--f  Vv-(OF)/dv 
me  J  me  meJ 


(10) 


where  Sv  is  the  boundary  of  the  velocity  space  at  infinity.  Since  no  particle  travels  at  an  infinite 
speed,  f  =  0  at  v  oo  and  the  contour  integral  in  (10)  vanishes.  As  a  result,  (10)  becomes 


[  0>  —  ■  Vvf  dv  =  -  —  [  Vv  ■  (cDF)  /  dv  =  -  —  (Vv  ■  (cDF)) 
J  me  me  J  me 


n 


n 


= - <Vvc D-F) - (d>  Vv  ■  F) 


me 
n  , 

= - (Vvd>  ■  F) 

mP 


mP 


(11) 


In  reaching  (11),  Vv  ■  F  =  0  is  applied  since  the  Lorentz  force  is  divergence-free  in  the  velocity 
space. 


The  RHS  of  (6)  is  denoted  as 

/  Hv)  (I) 


dv  = 


(  dn(<J>(v))\ 


V  dt  J  coll 

at  this  moment,  and  will  be  discussed  in  the  following  sections. 

Finally,  the  general  transport  equation  can  be  expressed  as 

dn(<J>)  n  /  dn(<J>)\ 

— ^  +  V  ■  (n(<J>v» - (F  ■  Vv<J>)  =  — ^ 

dt  me  \  at  J 


(12) 


(13) 


coll 


which  is  the  conservation  equation  for  the  density  of  the  macroscopically  averaged  quantity 

3.1.2  Five-Moment  Plasma  Fluid  Model 

Plasma  fluid  models  can  be  obtained  by  taking  moments  of  the  Boltzmann  equation.  In  this 
section,  the  five-moment  fluid  model  is  derived  by  taking  the  first  three  moments.  Higher-order 
fluid  models,  such  as  the  ten-moment  and  the  13-moment  models,  can  also  be  obtained  in  the  same 
manner. 

3. 1.2.1  The  Particle  Continuity  Equation 

In  (13),  setting  =  1  yields  the  particle  continuity  equation 

dn  .  .  fdn\ 


dn  „  .  / dn\ 

_+v.(„w)  =  y  =(v,-v> 


(14) 


coll 
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where  Vj  and  va  stand  for  the  ionization  and  attachment  frequencies  in  the  ionization  and 
attachment  collision  processes,  respectively,  which  can  be  physically  interpreted  as  the  number  of 
particles  created  and  annihilated  per  unit  time  per  unit  volume  due  to  collisions. 

3. 1.2.2  The  Momentum  Conservation  Equation 

Taking  d>  =  m  ev  yields  the  momentum  conservation  equation.  In  (13),  (<J>)  =  m  e(v) 
=  m  CU  and  (<f>v)  =  m  Q(vv).  Let  v  =  U  +  q  ,  where  q  stands  for  the  velocity  fluctuation 
around  the  mean  value  U,  we  have 


(®v)  =  me(vv)  =  me((UU)  +  ( Uq )  +  ( qU )  +  (qq))  =  me(UU  +  (qq)) 


(15) 


Defining  the  (symmetric)  pressure  tensor 

P  =  me  /  „  /  dv  =  me  f(v  -  V)(v  -  (/)  fir,  = 

yields 

1 

(d>v)  =  meUU  +  —  F 
n 

The  second  term  of  (13)  can  be  expressed  as 

V  ■  (n(d>v))  =  V  ■  (menUU  +  P)  =  meV  ■  (nUU  +  —  pj 


(16) 


(17) 


(18) 


In  the  third  term  of  (13),  Vvd>  =  m  eVvv  =  m  eI,  where  I  is  an  identity  tensor.  Therefore,  the 
third  term  becomes 


- (F  ■  Vv<J>)  = - ( F  ■  me\ I)  =  —n(qe(E  +  vxB ))  =  —qe(nE  +  nlJxB) 


mP 


mP 


(19) 


where  qe  denotes  the  electron  charge. 

Finally,  the  momentum  conservation  equation  is  obtained  as 

dnU  /  1  \  qe  1  fdnU\ 

— -  +  V-  (nUU  +  —  P)=—(nE  +  nUxB)+  —  (— — ) 
ot  \  me  )  me  me  \  ot  /con 

qP 

=  — ( nE  +  nUxB )  —  vcnU 
mP 


(20) 


where  vc  stands  for  the  total  collision  frequency. 
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3. 1.2.3  The  Energy  Conservation  Equation 


The  energy  conservation  equation  can  be  obtained  by  taking  d>  =  -mev2.  In  the  first  term  of 

(13),  (d>)  =  ^me(v2)  =  £  and  in  the  second  term,  n(d>v)  =  ^ men(v2v )  =  £  .  Using  the  definition 
of  the  average  quantity  (5),  we  have 


n(v2v)  =  J  v2v  f  dv  =  j(q  +  U)-(q  +  UXq  +  U)fdv 

=  [  q2q  f  dv  +  2U  ■  [  qq  f  dv  H - £U 

J  J  me 


(21) 


=  —  (Q  +  F-U  +  nU£ ) 


mP 


In  deriving  the  above  expression,  the  definition  of  the  pressure  tensor  (16)  has  been  employed,  and 
the  heat  flux  vector  is  defined  as 


Q  =  j  q2qfdv  =  ^  j \v  -  U\2(v  -  U)  f  dv 
As  a  result,  the  second  term  in  the  general  transport  equation  becomes 


(22) 


V  ■  (n(d>v»  =  V  ■  ( nU£  +  F-U  +  Q) 


(23) 


The  third  term  in  the  general  transport  equation  can  be  written  as 
n  , 

- (F  ■  Vvd>)  = 

me 

Finally,  the  energy  conservation  equation  can  be  expressed  as 

dn£  ,  / dn£\ 

— - 1-  V  ■  ( nU£  +  F  ■  U  +  Q)  =  qeE  ■  nU  +  ( — — )  =  qeE  ■  nU  —  Qen 

at  \  at  /co ii 


F-Vv(\mev2^  fdv  =  -jF-vfdv 
—  J  qe(E  +  vxB)  ■  v  f  dv  =  —  J  qeE  ■  v  f  dv 


—qeE  ■  nU 


(24) 


(25) 


where  Qe  stands  for  the  electron  energy  loss  frequency. 

Equations  (14),  (20),  and  (25)  are  the  governing  equations  of  the  five-moment  plasma  fluid 
model.  The  name  “five-moment”  refers  to  the  five  scalar  quantities  governed  by  the  model, 
which  are  n,  Ux,  Uy,  Uz,  and  £. 

It  should  be  noted  that  the  equations  obtained  from  the  general  transport  equation  are  not  closed. 
Due  to  the  second  term  V  ■  (n(<J>v)),  the  n-th  moment  equation  will  always  introduce  the  (n  +  1)- 
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th  macroscopic  moment.  As  a  result,  any  finite  set  of  moment  equations  have  more  unknowns 
than  equations.  To  obtain  a  closed  model,  some  additional  information,  limiting  assumption,  or 
additional  physical  setting,  is  always  needed. 

In  the  five-moment  model,  the  higher-order  moment,  the  heat  flux  vector  Q,  is  ignored,  and  the 
pressure  tensor  is  assumed  to  be  diagonal  and  isotropic  IP  =  PI,  where  P  is  the  scalar  pressure 
defined  as 


1  771  f  771  f 

P  =  -tr(lP)  =  — ^  tr[(v  —  U)(v  —  U)]  f  dv  =  — ^  \v  -  U\2  f  dv 
o  o  J  o  J 


(26) 


in  three  dimensions.  In  a  /V-dimensional  problem,  the  factor  1/3  becomes  1/iV.  In  the  above 
expression,  tr  stands  for  the  matrix  trace.  From  the  ideal  gas  law,  the  pressure  and  the 
temperature  are  related  by  the  equation  of  state  P  =  n  k  BT,  where  kB  =  8.617X10-5  eV/K  is 
the  Boltzmann  constant  and  T  is  the  temperature  in  Kelvin.  From  (26),  it  can  be  seen  that 


+  U2  —  v  ■  U  —  U  ■  v)  f  dv 


(27) 


Therefore,  the  pressure  and  the  energy  are  related  by 


3  1 

n£  =  -P  +  -menU2  (28) 

3  1  9 

Apparently,  n£  stands  for  the  total  energy  for  n  particles,  -P  and  -menl/z  stand  for  the  total 
internal  and  kinetic  energy  for  n  particles,  respectively.  In  general, 


n£  = 


P 

F7! 


+  -menU2 


(29) 


where  y  =  (iV  +  2)/iV  is  the  ratio  of  specific  heats,  which  equals  to  5/3  for  electrons  in  three 
dimensions  and  1 .4  for  air  particles. 


With  the  heat  flux  ignored  and  the  scalar  closure  between  the  pressure  and  the  energy  assumed, 
the  five-moment  plasma  fluid  model  is  finally  presented  as 


dnU 


dn 

—  +  V  ■  ( nU )  =  (l/j  -  va)n 


dt 


V-  (nUU  +  —  Pi)  =  -(nE  +  nUxB)-vcnU 
\  mP  /  mP 


dn£ 

dt 


+  V  ■  [U(n£  +  P)]  =  qeE  ■  nU  —  Qen 


(30) 

(31) 

(32) 
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3.1.3  Collision  Reactions,  Transport  Coefficients,  &  Electron  Energy  Distribution 
Function 

When  electrons  are  pushed  to  move  by  the  Lorentz  force,  they  collide  with  neutral  air 
particles  frequently.  Depending  on  the  energy  that  different  electrons  possess,  collision  events 
may  trigger  different  reactions.  Generally,  there  are  two  major  types  of  collisions,  namely,  elastic 
and  inelastic  collisions.  In  an  elastic  collision  event,  the  electron  collides  with  a  neutral  air 
particle  (a  nitrogen  or  an  oxygen  molecule,  for  example),  gets  bounced  to  some  scattering  angle 
depending  on  the  angle  of  incidence  and  losses  a  part  of  its  momentum  and  energy.  For  an 
inelastic  collision,  there  are  several  possibilities.  First,  the  electron  may  get  attached  to  (absorbed 
by)  the  air  particle  and  result  in  a  negative  ion,  which  is  called  an  attachment  reaction.  Second,  if 
the  electron  has  a  higher  energy,  it  can  excite  the  air  particle,  making  the  air  particle  rotate, 
vibrate,  or  get  excited  to  a  higher  energy  level.  Such  a  collision  process  is  called  an  excitation 
reaction.  Third,  if  the  electron  has  an  even  higher  energy,  it  is  able  to  kick  out  a  new  electron 
from  the  air  particle,  which  is  known  as  an  ionization  reaction.  In  this  work,  air  is  simply  treated 
as  20%  oxygen  plus  80%  nitrogen.  Various  collision  reactions  between  the  electron  and  the 
oxygen  and  nitrogen  molecules  are  summarized  in  Table  1  and  Table  2,  respectively.  In  these 
two  tables,  £k  stands  for  the  threshold  (minimum)  energy  for  the  k- th  collision  reaction  to  take 
place. 

Table  1.  Electron-Oxygen  Collision  Reactions  [31],  [32] 


Index 

Collision  Type 

Reaction  Formula 

(eV) 

01 

Three-body  attachment 

e  +  02  — »  (Tp;  O2*  +  02  — »  02  +  02 

0 

02 

Two-body  attachment 

6  +  O2  — *  0  +0 

0 

03 

Effective  momentum  transfer 

6  +  O2  — *  6  +  O2 

0 

04 

Rotational  excitation 

e  +  02  — *  e  +  02  (rot) 

0.02 

05 

Vibrational  excitation 

e  +  02  — >  e  +  02(v  =  1) 

0.19 

06 

Vibrational  excitation 

e  +  02  — *  e  +  02(v  =  Ires) 

0.19 

07 

Vibrational  excitation 

e  +  02  — »  e  +  02(v  =  2) 

0.38 

08 

Vibrational  excitation 

e  +  02  — *  e  +  02(v  =  2res) 

0.38 

09 

Vibrational  excitation 

e  +  02  — »  e  +  02(v  =  3) 

0.57 

010 

Vibrational  excitation 

e  +  02  — »  e  +  02(v  =  4) 

0.75 

Oil 

Metastable  excitation 

e  +  02  — *  e  +  02(a1Ag) 

0.977 

012 

Metastable  excitation 

e  +  02  — »  e  +  02(b1I]g) 

1.627 

013 

Metastable  excitation 

e  +  02  — >  e  +  02(c12;u  A+£u) 

4.5 

014 

Dissociative  excitation 

e  +  02  — »  e  +  0(3P)  +  0(3P) 

6.0 

015 

Dissociative  excitation 

e  +  02  — »  e  +  0(3P)  +  0(1D) 

8.4 

016 

Dissociative  excitation 

e  +  02  — »  e  +  0(1D)  +  0(1D) 

9.97 

017 

Ionization 

6  +  O2  — *  2e  +  Oj 

12.06 
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Table  2.  Electron-Nitrogen  Collision  Reactions  [33],  [34] 


Index 

Collision  Type 

Reaction  Formula 

^  (eV) 

N1 

Elastic  momentum  transfer 

e  +  N2  — »  e  +  N2 

0 

N2 

Rotational  excitation 

e  +  N2  — >  e  +  N2(rot) 

0.02 

N3 

Vibrational  excitation 

e  +  N2  — >  e  +  N2(v  =  Ires) 

0.29 

N4 

Vibrational  excitation 

e  +  N2  — »  e  +  N2(v  =  1) 

0.291 

N5 

Vibrational  excitation 

e  +  N2  — »  e  +  N2(v  =  2) 

0.59 

N6 

Vibrational  excitation 

e  +  N2  — »  e  +  N2(v  =  3) 

0.88 

N7 

Vibrational  excitation 

e  +  N2  — *  e  +  N2(v  =  4) 

1.17 

N8 

Vibrational  excitation 

e  +  N2  — >  e  +  N2(v  =  5) 

1.47 

N9 

Vibrational  excitation 

e  +  N2  — *  e  +  N2(v  =  6) 

1.76 

N10 

Vibrational  excitation 

e  +  N2  — »  e  +  N2(v  =  7) 

2.06 

Nil 

Vibrational  excitation 

e  +  N2  — >  e  +  N2(v  =  8) 

2.35 

N12 

Electronic  excitation 

e  +  N2  — »  e  +  N2(A3,  v  =  0  -  4) 

6.17 

N13 

Electronic  excitation 

e  +  N2  — »  e  +  N2(A3,  v  =  5  -  9) 

7.00 

N14 

Electronic  excitation 

e  +  N2  — »  e  +  N2(B3) 

7.35 

N15 

Electronic  excitation 

e  +  N2  — »  e  +  N2(W3) 

7.36 

N16 

Electronic  excitation 

e  +  N2  — »  e  +  N2(A3,  v  >  10) 

7.80 

N17 

Electronic  excitation 

e  +  N2  — »  e  +  N2(B'3) 

8.16 

N18 

Electronic  excitation 

e  +  N2  — »  e  +  N2(a'l) 

8.40 

N19 

Electronic  excitation 

e  +  N2  — >  e  +  N2(al) 

8.55 

N20 

Electronic  excitation 

e  +  N2  — *  e  +  N2(wl) 

8.89 

N21 

Electronic  excitation 

e  +  N2  — »  e  +  N2(C3) 

11.03 

N22 

Electronic  excitation 

e  +  N2  — »  e  +  N2(E3) 

11.87 

N23 

Electronic  excitation 

e  +  N2  — >  e  +  N2(a"l) 

12.25 

N24 

Electronic  excitation 

e  +  N2  — >  e  +  N2(sum  of  singlet  states) 

13.0 

N25 

Ionization 

e  +  N2  — »  2e  +  N2 

15.6 

The  various  collision  reactions  result  in  the  generation  and  annihilation  of  electrons  and  the 
momentum  and  energy  change  of  electrons,  which  are  accounted  for  by  the  transport  coefficients, 
including  the  total  ionization  coefficient  vb  the  total  attachment  coefficient  va,  the  total  collision 
frequency  vc,  and  the  energy  loss  frequency  Qe.  These  transport  coefficients  are  critical 
parameters  that  determine  the  behavior  of  the  breakdown  process  and  are  defined  as  the 
integration  of  the  corresponding  collision  cross  sections  multiplied  by  the  EEDF,  over  the  entire 
energy  spectrum  as  follows  [15]: 


9 


Approved  for  public  release;  distribution  is  unlimited. 


Vi  =  iVAir 


2 

77le 


sF0ds 


(33) 

(34) 

(35) 


Qe  ^Air 

N 

+^Air 


A  V  L  (W  +M^£\ 

me.  2-i  Mk  JEk  k[  k\  0  qe  del 


M 


k= elastic 

2 


d£ 


^  Z  £*(, 

/e=inelastic  k 


Xk(Tk£F0d£ 


(36) 


where  /VAir  denotes  the  number  density  of  the  ambient  air,  xfc  and  Mk  denote  the  mole  fraction  and 
the  mass  of  the  target  gas  species  of  collision  process  k,  respectively,  ak  denotes  the  cross  section 
of  collision  process  k,  and  F0  denotes  the  EEDF  normalized  by 

r 00  i 

I  £2F0(e)d£  =  1  (37) 

Jo 


_3 

which  has  units  of  eV  7. 

In  the  above  definitions,  the  collision  cross  sections  are  obtained  from  either  analytical 
calculations  or  experimental  measurements  and  are  usually  tabulated  as  functions  of  electron 
energy.  The  EEDF  is  usually  solved  from  the  Boltzmann  equation.  For  an  equilibrium  process, 
the  EEDF  has  a  Maxwellian  distribution  which,  when  normalized  by  (37),  becomes  a  straight 
line  in  a  logarithmic  plot,  with  the  slope  being  —l/kBT.  The  high-power  breakdown  process, 
however,  is  highly  non-equilibrium.  As  a  result,  the  Maxwellian  distribution  is  not  appropriate 
in  describing  the  electron  energy  distribution  in  such  a  process.  To  obtain  a  proper  distribution 
function,  the  classic  two-term  approximation  can  be  used  to  solve  the  Boltzmann  equation  taking 
into  consideration  the  breakdown  process,  where  the  electron  distribution  function  f  is  expanded 
in  terms  of  spherical  harmonics  as  [15] 


fir,  v,  t)  =  2-  & y  n(r,  t)  [F0  (f)  +  F1  (f)  cos  9] 


(38) 


where  F0  is  the  EEDF  representing  the  isotropic  part  of  the  electron  distribution  /  and  Ft  is  the 
anisotropic  perturbation  of  f.  Once  the  proper  EEDF  is  obtained,  it  can  be  used  to  integrate  with 
the  corresponding  cross  sections  to  calculate  the  transport  coefficients,  which  are  then  used  in  the 
five-moment  model  to  calculate  the  macroscopic  fluid  parameters. 
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3.2  Coupled  Electromagnetics-Plasma  System 

The  high-power  air  breakdown  process  can  be  described  by  the  following  coupled 
electromagnetic— plasma  system 


dnU 


dE 

e— - Vx/f  =  —qenll 

at 

dH 

/r  —  +  VxE  =  0 
at 
dn 

—  +  V  ■  ( nU )  =  Oi  -  va)n 


dt 


+  V  ■  (nUU  +  — PlA  =  —  (nE  +  nUxB)  —  vcnU 
\  mP  )  mP 


dn£ 

dt 


+  V  ■  [U(n£  +  P )]  =  qeE  ■  nU  —  Qen 


(39) 

(40) 

(41) 

(42) 

(43) 


In  the  above  equations,  e  and  q  denote  the  permittivity  and  permeability  of  the  medium, 
respectively.  The  electric  and  magnetic  fields  E  and  H  are  governed  by  Maxwell’s  equations  (39) 
and  (40),  while  the  electron  density  n,  mean  velocity  U,  mean  energy  £  and  pressure  P  are 
governed  by  the  five-moment  plasma  fluid  equations  (41)— (43),  which  are  essentially  the 
compressible  Euler  equations  in  fluid  dynamics.  As  described  in  the  preceding  section,  in  order 
to  close  the  system,  the  heat  flux  is  neglected  and  the  scalar  closure  is  applied  to  relate  the  energy 

p  1  o 

and  the  pressure  as  n£  =  j-  +  -menU  . 


3.3  The  DGTD  Solution  of  the  Coupled  System 

To  solve  the  coupled  equations  using  the  DGTD  method,  both  Maxwell’s  and  the  Euler 
equations  are  first  written  uniformly  into  the  following  conservation  form 

dG 

+  V-g(G)=5(G)  (44) 


where  D  denotes  the  material  parameter  tensor,  G  denotes  the  unknown  quantity  vector  which 
consists  of  the  conservative  variables,  $  denotes  the  physical  flux  tensor  and  5  denotes  the 
source  term. 


The  strong  form  of  these  governing  equations  can  then  be  expressed  as 

h  (8*  ~  5)  ■  n  dS  +  f  lt  5(G)  dE  (45) 

Jve 

where  stands  for  the  testing  function  defined  in  the  tetrahedral  element  V e  and  ($*  —  5  )  1  n  ~ 
refers  to  the  total  flux  defined  on  the  boundary  of  each  tetrahedral  element. 


11 


Approved  for  public  release;  distribution  is  unlimited. 


3.3.1  Maxwell’s  Equations 

For  Maxwell’s  equations,  the  material  parameter  tensor  is 

D  =  diag{e,  e,  e,  /t,  /t,  jU>  (46) 

the  conservative  variables  are  simply  chosen  as  (with  T  being  the  transpose  operator) 

G=(E1,E2>E3,H1,H2>H3y  (47) 

the  flux  tensor  is 


0 

-h3 

H  2 

H3 

0 

-h2 

0 

0 

E3 

E2 

-Es 

0 

Ei 

e2 

~Et 

0 

and  the  source  term  is 

S  =  [—qenUx,  —qenUy,  —qenUz,  0,0,0)T  (49) 

Due  to  their  hyperbolic  property,  the  upwind  flux  can  be  used  for  Maxwell’s  equations  [14],  [18] 


(r-S)-n 


^y(r+  nx[£]+nxnx[/f]) 

-^(Z+nxm-nxnxm) 


In  the  above  expression,  the  addition  (■)  and  jump  [■]  functions  are  defined  as  follows 

(a)  =  a+  +  a~,  [a]  =  a+  —  a~  (51) 

where  the  superscripts  —  and  +  indicate  the  inside  and  the  outside  of  a  surface,  respectively. 
Note  that  these  definitions  apply  to  both  scalar  and  vector  variables. 

3.3.2  Euler’s  Equations 

For  Euler’s  equations,  X)  is  an  identity  tensor,  the  conservative  variables  are  chosen  as 

G  =  (n,nUx,nUy,  nUz,  nE)  (52) 
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Once  these  conservative  variables  are  obtained,  the  primitive  variables  can  be  obtained  as  U  = 
nU/n,  £  =  n£/n,  and  P  =  (y  —  1)  (n£  —  ^menU2^j,  which  are  needed  in  the  calculation  of  the 
flux  tensor 
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(53) 


and  the  source  term 


[Vi(£)  -  va(£)]n 


Re 


S  =  —  [n(E  +  Fnc)  +  nUx(B  +  Bmc )]  -  vc(£)nU 

me 

qe(E  +  Fnc)  ■  nU  -  Qe(£)n 


(54) 


Here,  £'inc  and  Bmc  stand  for  the  electric  field  and  magnetic  flux  of  the  incident  high-power 
pulse,  respectively.  The  transport  coefficients  are  all  functions  of  the  mean  energy  £  and  can  be 
obtained  by  interpolating  the  tabulated  values  calculated  using  the  method  introduced  in  the 
preceding  section. 


To  solve  the  strong  form  (45)  of  the  Euler  equations,  the  local  Lax-Friedrichs  flux  [18] 


(r-S)-n  =  -ra-n-A[G]} 


can  be  used,  where 


A  =  max 

ge{G+,G-} 


jl</(0)l  + 


(55) 


(56) 


is  the  characteristic  velocity  of  Euler’s  equations.  The  resulting  ordinary  differential  equations 
from  both  Maxwell’s  and  Euler’s  equations  are  integrated  in  time  using  the  Runge-Kutta  method 
to  obtain  the  time-domain  responses  and  the  coupling  between  these  two  physics  is  carried  out 
through  the  source  terms  (49)  and  (54). 
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4.0  RESULTS  AND  DISCUSSION 


In  this  section,  several  examples  are  given  to  illustrate  the  EEDF  and  transport  coefficients 
and  to  demonstrate  the  breakdown  process  and  behaviors  under  different  conditions.  In  all  the 
numerical  examples,  the  gas  under  consideration  is  ambient  air,  which  is  assumed  to  consist  of 
80%  nitrogen  and  20%  oxygen,  with  in  total  42  different  (elastic  and  inelastic)  collision  processes 
listed  in  Table  1  and  Table  2. 

4.1  EEDF  &  Transport  Coefficients 

First,  the  Maxwellian  and  non-Maxwellian  EEDFs  at  different  mean  energy  are  shown  in  this 
section  to  demonstrate  their  distinct  difference.  As  seen  in  Figure  1,  the  EEDFs  with  the 
Maxwellian  distribution  are  all  straight  lines  in  the  semi-log  plot,  indicating  that  the  number  of 
electrons  decreases  exponentially  as  the  energy  increases.  The  non-Maxwellian  EEDFs  obtained 
by  solving  the  Boltzmann  equation  using  the  BOLSIG+  package  [15]  however,  decrease  much 
faster  as  the  energy  increases,  which  indicates  that  fewer  electrons  have  higher  energies  in  a  non¬ 
equilibrium  case.  The  higher  energy  tail  of  the  Maxwellian  EEDF  results  in  a  larger  total 
ionization  frequency  v{  —  v  a  when  the  mean  energy  is  higher  than  1  eV,  as  shown  in  Figure  2a. 
The  total  collision  (Figure  2b)  and  power  loss  (Figure  2c)  frequencies  for  the  Maxwellian  EEDF 
are  also  larger  than  the  non-Maxwellian  EEDF  at  lower  mean  energies. 


(a)  (b) 

Figure  1.  The  Maxwellian  and  Non-Maxwellian  EEDFs  at  a  Mean  Energy  of  (a)  1  eV  and 

(b)  10  eV 
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Figure  2.  The  Transport  Coefficients  as  Functions  of  the  Mean  Energy,  Calculated 
Using  the  Maxwellian  and  Non-Maxwellian  EEDFs.  (a)  Total  Ionization  Frequency,  (b) 
Total  Collision  Frequency  and  (c)  Power  Loss  Frequency 


4.2  Air  Breakdown  Using  Different  EEDFs 

To  demonstrate  the  effects  of  different  EEDFs  in  the  prediction  of  air  breakdown,  a  numerical 
example  is  presented  in  this  section,  where  a  set  of  5. 5 -MV/m,  200-MHz  and  y-polarized 
modulated  difference-of-double-exponentials  (DEXP)  pulse  [37]  traveling  through  air  is 
considered.  As  illustrated  in  Figure  3,  the  simulation  domain  is  30  m  long  in  the  z-direction 
which  is  truncated  with  the  absorbing  boundary  condition  (ABC)  and  0.5  m  in  both  the  x-  and  the 
y-directions  which  are  truncated  with  the  perfect  electric  and  magnetic  conductors,  respectively. 
The  initial  electron  density  in  the  air  is  set  as  106/m3.  The  ambient  temperature  of  the  air  is  set 
as  the  room  temperature  at  300  K  and  the  ambient  pressures  under  consideration  is  760  Torr. 
Figure  4  presents  the  simulated  results  using  the  Maxwellian  and  non-Maxwellian  EEDFs. 
Clearly,  the  electric  field  distribution  obtained  using  the  non-Maxwellian  EEDF  (Figure  4a) 
remains  undisturbed,  which  is  the  same  as  the  pulse  traveling  in  the  vacuum.  However,  when  the 
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Maxwellian  EEDF  is  used  to  estimate  the  transport  coefficients,  the  higher  energy  tail  results  in 
larger  ionization  coefficient  which  leads  to  a  larger  increase  of  the  electron  density,  as  shown  in 
Figure  4c.  The  electron  density  obtained  by  the  Maxwellian  EEDF  on  the  plateau  is  about  100 
times  larger  than  that  obtained  by  the  non-Maxwellian  EEDF.  This  results  in  a  much  larger 
absorption  of  the  incident  field  and  leads  to  the  tail  erosion  of  the  incident  HPM  pulse,  as  can  be 
seen  in  Figure  4a.  Other  physical  parameters,  including  the  electron  velocity,  energy  and 
temperature,  also  differ  significantly  between  the  results  obtained  from  the  Maxwellian  and  non- 
Maxwellian  distributions.  Shown  in  Figure  4b  and  Figure  4e  are  the  z-component  (the 
longitudinal  component)  of  the  electric  field  and  electron  velocity,  respectively.  This  component 
is  present  because  of  the  magnetic  Lorentz  force  term  qenUxB.  However,  since  the  Uy 
component  is  very  small  compared  to  the  speed  of  light,  the  longitudinal  component  is  very 
small.  But  it  will  become  significant  when  the  kinetic  energy  of  the  electrons  is  high  enough 
when  compared  to  the  speed  of  light.  Moreover,  it  is  interesting  to  see  that  both  the  electric  field 
Ez  and  the  electron  velocity  Uz  are  positive,  which  is  distinctively  different  from  Ey  and  Uy. 
This  is  due  to  the  fact  that  the  Uy  and  Bx  components  are  always  in-phase.  Resulted  from  their 
cross  product,  the  Uz  component  becomes  purely  positive  and  so  does  Ez.  The  positive  motion 
pushes  the  electrons  toward  the  front  edge  of  the  plasma  and  (partially)  results  in  the  large 
density  gradient  on  the  plasma  edge.  Also,  it  is  very  clear  that  the  oscillating  frequency  of  Uz  is 
twice  as  Uy ,  due  to  the  higher-order  oscillation  generated  by  the  magnetic  Lorentz  force. 


Figure  3.  Illustration  of  the  Solution  Domain 
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Figure  4.  The  Electric  Field,  Electron  Density,  Velocity,  Total  Energy  and  Temperature 
Distributions  in  Space  as  the  Modulated  DEXP  Pulse  Travels  Through  the  300-K,  760-Torr 
Air.  Comparisons  Are  Made  Between  Using  the  Maxwellian  and  Non-Maxwellian  EEDFs 
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4.3  Air  Breakdown  Versus  Ambient  Pressure 


To  investigate  the  breakdown  behaviors  under  different  ambient  pressures,  a  set  of  300-kV/m, 
200-MHz  modulated  DEXP  pulse  traveling  through  air  is  considered.  As  in  the  previous 
example,  the  initial  electron  density  in  the  air  is  set  as  106/m3  and  the  ambient  temperature  of 
the  air  is  set  as  the  room  temperature  at  300  K.  The  ambient  pressures  under  consideration  are 
600  Torr  and  10  Torr,  which  are  referred  to  as  the  high-  and  the  low-pressure  cases, 
respectively.  Figure  5  and  Figure  6  present  the  simulated  results  for  these  two  cases, 
respectively.  In  the  high-pressure  case,  the  air  does  not  break  down  (see  Figure  5c)  and  the  HPM 
pulse  travels  through  the  air  undisturbed.  This  can  be  seen  clearly  in  Figure  5a  where  the  electric 
field  distribution  in  the  air  is  the  same  as  that  in  the  vacuum  where  no  breakdown  can  take  place. 
In  the  low-pressure  case,  however,  air  breakdown  takes  place  and  the  tail  of  the  HPM  pulse 
cannot  travel  through  the  air,  as  can  be  seen  in  Figure  6a  and  Figure  6c.  The  effect  of  pressure  on 
the  air  breakdown  process  can  be  explained  as  follows.  In  the  high-pressure  case,  due  to  the 
much  higher  air  particle  density,  the  electrons  collide  much  more  frequently  with  air  particles, 
which  causes  more  electron  momentum  and  energy  loss.  Therefore,  it  is  more  difficult  for  the 
electron  density  to  increase  in  the  high-pressure  case  than  in  the  low-pressure  case.  Compared  to 
the  high  pressure  results  shown  in  Figure  5d  and  Figure  5f,  the  electron  velocity  and  energy 
shown  in  Figure  6d  and  Figure  6f  are  increased  to  a  much  higher  level,  which  leads  to  a  much 
higher  ionization  frequency  that  causes  the  breakdown.  This  example  has  been  reported  in  [38], 
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Figure  5.  The  Electric  Field,  Electron  Density,  Velocity,  Total  Energy  and  Temperature 
Distributions  in  Space  as  the  Modulated  DEXP  Pulse  Travels  Through  the  300-K,  600-Torr 

Air 
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Figure  6.  The  Electric  Field,  Electron  Density,  Velocity,  Total  Energy  and  Temperature 
Distributions  in  Space  as  the  Modulated  DEXP  Pulse  Travels  Through  the  300-K,  10-Torr 

Air 
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Next,  the  breakdown  time  as  a  function  of  ambient  pressure  is  calculated  and  presented  in 
Figure  7.  The  breakdown  time  is  usually  defined  as  the  time  needed  for  the  electron  density  to 
increase  to  108  times  its  initial  value  [39],  In  this  investigation,  a  4.23-MV/m,  2.82-GHz 
sinusoidal  wave  travels  through  air  and  causes  the  breakdown.  The  initial  electron  density  in  the 
air  is  set  as  106/m3  and  the  ambient  temperature  of  the  air  is  set  as  the  room  temperature  at  300 
K.  In  Figure  7,  the  air  breakdown  time  versus  ambient  pressure  is  presented  and  compared  with 
reference  results  obtained  by  a  particle-in-cell  with  Monte-Carlo  collision  (PIC-MCC)  simulation 
[39],  [40]  and  a  fluid  simulation  with  a  modified  EEDF  [9].  Excellent  agreement  is  observed, 
which  validates  our  implementation  of  the  coupled  DGTD  solver. 


Pressure  (Torr) 

Figure  7.  Breakdown  Time  in  Air  at  4.23  MV/m  and  2.82  GHz.  The  Results  From  PIC- 
MCC  Simulations  and  a  Fluid  Model  With  a  Modified  EEDF  Are  Provided  as  References 


4.4  Plasma  Shielding  Effect 

Air  breakdown  and  plasma  shielding  triggered  by  geometry  is  demonstrated  in  this  example. 
As  shown  in  Figure  8a,  the  solution  domain  considered  in  this  example  is  a  parallel  plate 
waveguide  with  a  metallic  wall  placed  between  the  two  parallel  plates,  which  forms  a  rectangular 
aperture  (denoted  in  red  in  the  figure).  The  solution  domain  is  truncated  from  the  left  and  the 
right  using  the  ABC.  A  25-GHz,  2.0-MV/m  and  vertically  (the  y  direction)  polarized  plane  wave 
with  a  tapered  sinusoidal  temporal  profile  is  launched  from  the  left  boundary  and  propagates 
toward  the  right  direction  (the  z  direction).  A  100-torr,  300-K  air  is  assumed  to  be  confined  in 
the  aperture  area,  with  the  rest  of  the  domain  filled  with  vacuum.  The  initial  spatial  profile  of  the 
electron  density  is  set  to  be  sinusoidal  in  both  y  and  z  directions,  which  reaches  its  maximum 
value  of  1015/m3  at  the  center  of  the  aperture.  Two  observation  points  PI  {x  =  y  =  0 . 0  mm, 
z  =  0.0  mm)  and  P2  (x  =  y  =  0 . 0  mm,  z  =  2.81  mm)  are  also  shown  in  Figure  8a. 

Besides  the  five-moment  model  presented  in  this  report,  this  example  has  been  simulated  using 
a  simplified  diffusion  model  [14],  where  only  the  electron  density  and  velocity  were  considered. 
It  should  be  pointed  out  that  in  the  diffusion  model,  the  plasma  parameters  (ionization,  attachment, 
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and  momentum  transfer  frequencies,  etc.)  were  obtained  by  empirical  models  based  on  measured 
data,  which  were  nonlinear  functions  of  the  reduced  effective  electric  field  Ee{{/P,  whereas  in  the 
five-moment  model,  the  plasma  parameters  are  obtained  by  integrating  collision  cross  sections 
with  non-Maxwellian  EEDF,  which  are  nonlinear  functions  of  the  electron  energy.  As  a  result, 
these  two  models  are  very  different  in  terms  of  the  physical  quantities  under  consideration, 
governing  equations  and  plasma  parameters  calculation. 

Shown  in  Figure  8b,  Figure  8c  and  Figure  8e  are  the  electric  fields  recorded  at  PI  and  P2  and 
the  electron  density  evolution  recorded  at  PI.  The  comparisons  are  made  between  these  two 
different  models.  Excellent  agreements  are  achieved,  which  show  that  the  physical  quantities  of 
interest  can  be  obtained  by  very  different  physical  models.  Shown  in  Figure  8d  is  the  comparison 
of  electric  fields  recorded  at  PI  and  P2.  Apparently,  after  breakdown  takes  place  at  0.5  ns,  the 
incident  fields  are  blocked  and  only  a  small  portion  of  energy  can  transmit  to  the  observation 
point  P2.  In  Figure  8f,  the  electron  energy  and  temperature  at  PI  with  respect  to  time  are  given, 
which  provide  more  physical  insight  of  the  process  and  fidelity  of  the  simulation. 
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Figure  8.  Microwave  Breakdown  and  Plasma  Shielding  in  a  Metallic  Aperture,  (a) 
Illustration  of  the  Solution  Domain,  (b)  and  (c)  Comparison  of  the  Electric  Fields 
Simulated  Using  the  Diffusion  Model  and  the  Five-Moment  Model,  Recorded  at  PI  and 
P2,  Respectively,  (d)  Temporal  Response  of  the  Electric  Fields  Recorded  at  PI  and  P2.  (e) 
Temporal  Evolution  of  the  Electron  Density  Recorded  at  PI.  (f)  Temporal  Evolution  of 
the  Electron  Energy  and  Temperature  Recorded  at  PI 
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4.5  Air  Breakdown  Around  PEC  Objects 


As  the  last  example,  air  breakdown  around  two  perfectly  electric  conducting  (PEC)  cylinders 
are  simulated  using  the  five-moment  model.  This  example  is  designed  to  demonstrate  the  plasma 
pattern  formation  and  the  scattering  characteristic  change  after  air  breakdown.  Here,  two  PEC 
cylinders  with  the  same  radius  of  0.25  mm  are  placed  in  free  space.  A  5.5-MV/m,  200-GHz  and 
y-polarized  plane  wave  is  incident  from  the  z-direction.  Shown  in  Figure  9  is  the  scattered 
electric  field  distribution  when  the  entire  simulation  domain  is  filled  with  vacuum,  which  serves 
as  the  baseline  for  comparison.  When  the  space  is  filled  with  100-Torr,  300-K  ambient  air  with 
106/m3  initial  electrons,  the  incident  HPM  triggers  breakdown,  first  between  the  two  PEC 
cylinders  where  the  field  is  intensified.  Due  to  the  mutual  couplings  among  the 
electromagnetic  fields,  the  PEC  scatterers  and  the  plasma  fluid,  the  ionized  electrons 
follow  a  similar  distribution  as  the  scattered  field  pattern,  as  shown  in  Figure  10c,  Figure  lOe  and 
Figure  lOf.  At  0.25  ns,  although  increased  by  1 1  orders  of  magnitude,  the  electron  density  is  still 
low,  with  a  maximum  density  being  1.23xl017/m3.  The  electron  oscillation  does  not  contribute 
significantly  to  the  scattered  electric  field.  As  a  result,  the  electric  fields  presented  in  Figure  10a 
and  Figure  10b  are  very  similar  to  those  given  in  Figure  9.  As  the  breakdown  process  continues, 
the  electron  density  and  energy  increase  further.  At  0.5  ns,  the  maximum  electron  density  has 
increased  by  another  5  orders  of  magnitude,  reaching  as  high  as  5.92xl022/m3  (Figure  11c). 
The  strong  electron  oscillation  starts  to  alter  the  scattered  field  pattern.  As  can  be  seen  in  Figure 
11a  and  Figure  lib,  the  scattered  fields  are  significantly  different  from  those  recorded  at  0.25  ns 
and  those  in  the  vacuum  case.  Due  to  the  strong  electron  radiation,  certain  micro-patterns  are 
formed  in  Figure  11,  which  have  a  characteristic  length  much  smaller  than  the  wavelength  of  the 
incident  electromagnetic  field.  This  is  an  illustration  of  the  multiscale  feature  caused  by  the 
complicated  electromagnetic— plasma  interactions,  which  can  only  be  captured  by  a  numerical 
method  with  a  high  spatial  resolution. 
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Scattered  Electric  Field  Distribution  in  Vacuum,  (a)  Ey;  (b)  Ez 
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Time:  0.25  ns 
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Figure  10.  Electric  Field  and  Plasma  Quantities  Observed  at  0.25  ns.  (a)  Ey;  (b)  Ez ;  (c) 
Electron  Density;  (d)  The  Magnitude  of  the  Electron  Velocity;  (e)  Electron  Energy  and 

(f)  Electron  Temperature 
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Figure  11.  Electric  Field  and  Plasma  Quantities  Observed  at  0.5  ns.  (a)  Ey ;  (b)  Ez;  (c) 
Electron  Density;  (d)  The  Magnitude  of  the  Electron  Velocity;  (e)  Electron  Energy  and 

(f)  Electron  Temperature 


To  better  understand  the  physical  process,  the  plasma  quantities  and  secondary  fields  are  plotted 
along  the  z-axis  in  Figure  12.  As  shown  in  Figure  12a,  the  electron  density  follows  the  scattered 
wave  pattern  in  the  lit  region  (z  <  0  mm).  At  0.25  ns,  the  electron  energy  has  a  highest  peak  at 
around  z  =  0  mm  (Figure  12b),  which  leads  to  the  fastest  increase  of  the  electron  density  in  this 
region.  When  its  density  reaches  1022/m3  at  0.5  ns,  the  electrons  start  to  shield  the 
electromagnetic  field.  This  results  in  a  strong  reflection  of  the  field  and  leads  to  a  higher  energy 
peak  at  around  z  =  —  0.  6mm,  which  will  cause  a  stronger  ionization  and  faster  increase  of  the 
electron  density  in  that  region.  When  this  process  continues,  density  spikes  will  gradually  form 
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and  propagate  toward  the  source  of  the  incident  field.  This  phenomenon  is  presented  in  Figure 
13  and  Figure  14  in  two  and  one  dimension,  respectively,  where  the  process  of  the  plasma  pattern 
formation  can  be  observed  very  clearly.  Such  a  filamentary  plasma  array  formation  was  also 
observed  in  experiments  [3]-[5],  in  which  a  similar  plasma  pattern  has  been  recorded  by  a  high¬ 
speed  camera.  From  these  two  figures,  the  propagation  velocity  of  the  plasma  filaments  can  be 
estimated  to  be  about  1500  km/s,  which  is  on  the  same  order  of  the  theoretic  estimation  of  2yfD~vi 
[6]  where  De  and  v,  are  the  free  electron  diffusion  and  the  ionization  coefficients,  respectively. 
From  Figure  15,  the  spacing  between  two  adjacent  filaments  can  be  seen  very  clearly,  which  is 
about  a  quarter  wavelength  A0/4  =  0.375  mm  in  free  space.  This  is  because  the  strong  field 
reflection  from  the  filaments  results  in  a  standing  wave  pattern  in  the  lit  region,  which  generates 
plasma  filaments  with  a  similar  standing  wave  pattern. 
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Ey  (MV/m)  Temperature  (eV)  Density  (m*3) 


z-axis  (mm) 


(a) 
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z-axis  (mm) 


z-axis  (mm) 

(b) 


(d) 


(e)  (f) 

Figure  12.  Comparisons  of  Plasma  Quantities  and  the  Secondary  Electromagnetic  Fields 
Observed  at  0.25  and  0.50  ns  Along  the  z-Axis.  (a)  Electron  Density;  (b)  Electron  Energy; 
(c)  Electron  Temperature;  (d)  Electron  Velocity  Uy ;  (e)  Electric  Field  Ey  and  (f)  Magnetic 

Field  Hx 
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Time:  0.25  ns 


Time:  0.50  ns 


Time:  0.75  ns 


Time:  1.00  ns 


(e)  (f)  (g)  (h) 

Figure  13.  Electron  Density  Distribution  in  Linear  Scale,  Observed  at  (a)  0.25  ns;  (b)  0.50 
ns;  (c)  0.75  ns;  (d)  1.00  ns;  (e)  1.25  ns;  (f)  1.50  ns;  (g)  1.75  ns;  and  h)  2.00  ns 


29 


Approved  for  public  release;  distribution  is  unlimited. 


( (((13 


Density  (xIO22  m'3) 


I  1  I  1  I  1  I 


Time:  0. 75  ns| 


Time:  1.25  ns| 

A/A- 


I  '  I  '  I  1  I 


I  1  I  1  I  1  I 


1  I  1  I  1  1  1 


I  1  I  '  I  1  I 


7/me:  1.50  ns| 

aaA-,1 


'  I  '  I  '  I  '  I  ' 


7/me:  1.75  ns| 


— -AAAA 


I  '  I  '  I  '  I 


±  *  ■  «  ■— t 


7/me:  2.00  ns| 


LAAAAA 


-2.50  -1.25  0.00  1.25  2.50  -2.50  -1.25  0.00  1.25  2.50 

z-axis  (mm)  z-axis  (mm) 

Figure  14.  Electron  Density  Distribution  Along  the  z-Axis,  in  Linear  Scale,  Observed  at 

Different  Time  Instants 


(a)  (b) 

Figure  15.  Spacing  of  the  Plasma  Filaments,  (a)  2D  Plot;  (b)  ID  Plot 
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It  is  shown  in  this  section  that  this  very  interesting  and  highly  complicated  physical 
phenomenon  can  be  reproduced  successfully  with  the  proposed  numerical  method,  through  which 
many  physical  details  can  be  revealed. 

5.0  CONCLUSIONS 

In  this  report,  the  five -moment  plasma  fluid  model  has  been  employed  to  develop  a  DGTD- 
based  simulation  of  the  HPM  breakdown  in  air,  which  is  able  to  provide  more  physical  insight 
with  a  higher  fidelity  compared  to  simpler  models.  To  characterize  of  the  highly  non-equilibrium 
breakdown  process,  the  transport  coefficients  are  calculated  using  the  non-Maxwe Ilian  EEDFs 
solved  from  the  BOLSIG+  package.  The  self-consistent  Maxwell— Euler  system  equations  have 
been  solved  using  a  coupled  DGTD  method,  where  the  upwind  flux  is  used  for  Maxwell’s 
equations  and  the  Lax-Friedrichs  flux  is  used  for  Euler’s  equations.  The  employment  of  the  DGTD 
method  with  higher-order  nodal  basis  functions  provides  a  sufficient  spatial  resolution  to  capture 
the  fast-varying  micro-patterns  of  the  electromagnetic  fields  and  the  plasma  fluids.  Several 
numerical  examples  have  been  presented  to  show  the  difference  between  the  use  of  the  Maxwellian 
and  non-Maxwellian  EEDFs,  illustrate  the  air  breakdown  behavior  under  different  ambient 
pressures,  and  demonstrate  the  plasma  formation  and  electromagnetic  shielding  phenomena.  With 
the  method  presented  in  this  report,  many  interesting  physical  phenomena  can  be  modeled  and 
simulated  with  a  high  fidelity. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS  AND  ACRONYMS 
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De 

DEXP 

DGTD 
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EEDF 

Eeff/P 

Emc 
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Fo 

FDTD 

FETD 
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GPU 

1/ 

HPM 
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kB 

LHS 

Mk 

me 

^Air 
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PEC 

PIC-MCC 


<2 

<2e 


Qe 

RHS 


r  =  (x,y,z) 

S 

T 

U 

K 

V  Vy,  VZ ^ 


magnetic  flux  of  the  incident  high-power  pulse 

material  parameter  tensor 

free  electron  diffusion  coefficient 

difference  of  double  exponentials 

discontinuous  Galerkin  time-domain 

mean  energy 

electric  field 

electron  energy  distribution  function 

reduced  effective  electric  field 

electric  field  of  the  incident  high-power  pulse 

physical  flux  tensor 

Lorentz  force 

EEDF 

finite-difference  time-domain 

finite-element  time-domain 

electron  distribution  function 

conservative  variables  (unknown  quantity  vector) 

graphic  processing  unit 

magnetic  field 

high-power  microwave 

identity  tensor 

Boltzmann  constant 

left-hand  side 

testing  function 

mass  of  the  target  gas  species  of  collision  process  k 

electron  mass  at  rest 

number  density  of  the  ambient  air 

particle  density 

pressure  tensor 

scalar  pressure 

perfectly  electric  conducting 

particle-in-cell  with  Monte-Carlo  collision 

heat  flux  vector 

electron  energy  loss  frequency 

electron  charge 

right-hand  side 

physical  space 

source  term  of  conservation  equations 

temperature 

mean  velocity 

tetrahedral  element 

velocity  space 

mole  fraction  of  the  target  gas  species  of  collision  process  k 
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£k 
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Aq 
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va 

vc 

ak 
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threshold  energy 
permittivity  of  the  medium 
macroscopically  averaged  quantity 
ratio  of  specific  heats 
characteristic  velocity 
wavelength  in  free  space 
permeability  of  the  medium 
attachment  frequency 
total  collision  frequency 
ionization  frequency 
cross  section  of  collision  process  k 
total  flux 
addition  function 
jump  function 
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